
clear;
close all;
%% select channel
P="D:\APPdocu\python\UWA_GAN\data";
GAN_Path=P;
flag=1; % 1: NOF1; 2: NSC1
channel1='Generated_data\NOF1_256';
channel2='Generated_data\NCS1_256'; 
filepath1 = fullfile(GAN_Path, channel1, 'mat');
filepath2 = fullfile(GAN_Path, channel2, 'mat');
% 
% if flag==1
%    fs_tau=16000;
%    fs_t=fs_tau/510/4;
% else
%    fs_tau=16000;
%    fs_t=fs_tau/2048;
% end

temp=22;
filename1 = fullfile(filepath1, ['mat_' sprintf('%03d', temp) '.mat']);
filename2 = fullfile(filepath2, ['mat_' sprintf('%03d', temp) '.mat']);

var_struct1=load(filename1);
name_cell=fieldnames(var_struct1);
mat1=getfield(var_struct1,char(name_cell));


var_struct1=load(filename2);
name_cell=fieldnames(var_struct1);
mat2=getfield(var_struct1,char(name_cell));


h1_real=mat1(:,:,1);
h1_img=mat1(:,:,2);
h1=h1_real+1i*h1_img;

h2_real=mat2(:,:,1);
h2_img=mat2(:,:,2);
h2=h2_real+1i*h2_img;


%% 定义信道参数
%--------------------------------------------------------------------------

M_mod = 4;                         % 4QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 
fc = 12500;
c = 1500;                                                     % speed of light (m/s)
N = ofdmSym;

Fn = dftmtx(N);       % Generate the DFT matrix
Fn = Fn./norm(Fn);    % normalize the DFT matrix

Fm = dftmtx(M);       % Generate the DFT matrix
Fm = Fm./norm(Fm);    % normalize the DFT matrix

Md = M-Mzp;

maxDoppler=20;       % Maximum normalised digital velocity
c1_AFDM = (2*(maxDoppler) + 1)/(2*M); % Satisfying the orthogonality condition
c2_AFDM = 1/(2*M); 

%% 时延降采样且时间升采样
[~,G1]=channel_acc(h1,M,N,1);
[~,G2]=channel_acc(h2,M,N,2);
%% 发送信号
data_grid = zeros(M,N);
data_grid(1:Md,:) = 1;
N_syms_perfram = sum(sum(data_grid));

delta_f = scs;
T = 1/delta_f;
             
%%                       OTFS BER Calculation

     data_na= randi([0,1],N_syms_perfram*k,1);
     data=qammod(reshape(data_na,k,Md*N), M_mod,'gray','InputType','bit');
     X = Gen_2D_data_grid(N,M,data,data_grid);
    
    
    %% OTFS modulation%%%%
        X1_tilda=AFDMmod(X,c1_AFDM,c2_AFDM);
        s1= reshape(X1_tilda,N*M,1);
    
     %% channel output%%%%%

    r1=G1*s1;
    r2=G2*s1;
 
     %% OTFS demodulation%%%%
    Y1_tilda=reshape(r1,M,N);
    Y1 = Y1_tilda*Fn;
    y1=reshape(Y1.',N*M,1);
%% detect
   eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
    SNR_dB = -6:3:15;
    SNR = 10.^(SNR_dB/10);
    noise_var_sqrt = sqrt(1./SNR);
    sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;     

    errorRate1 = comm.ErrorRate('ResetInputPort',true);
    errorRate2 = comm.ErrorRate('ResetInputPort',true);

    berOTFS1 = zeros(length(SNR),3); 
    berOTFS2 = zeros(length(SNR),3);

    for m = 1:length(SNR)
        disp(m)
        noise= sqrt(sigma_2(m)/2)*(randn(size(s1)) + 1i*randn(size(s1)));
        r1_n=r1+noise;
        r2_n=r2+noise;

        [est_LMMSE_r1_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r1_n,G1,c1_AFDM,c2_AFDM); 
        [est_LMMSE_r2_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r2_n,G2,c1_AFDM,c2_AFDM);
        berOTFS1(m,:)=errorRate1(data_na,est_LMMSE_r1_n,1).';
        berOTFS2(m,:)=errorRate1(data_na,est_LMMSE_r2_n,1).';
    end

%% figure

figure

semilogy(SNR_dB,berOTFS1(:,1),'Marker','o','LineWidth',2,'Color',[0 0 1]);             %Plot simulated BER w/ OTFS
hold on;
semilogy(SNR_dB,berOTFS2(:,1),'Marker','*','LineWidth',2,'Color',[1 0 0]);
% semilogy(SNR,berMChannel(:,1),'--b');             %Plot simulated BER w/ C-OTFS
    

ylabel('BER');
xlabel('SNR/dB');
legend( 'NOF1',"NCS1");
title("AFDM输出误码率")
xlim([-6,15]);
grid on;
hold off;
